function C = D1(C,CON)
%*******************************************************************
%*** * Equations of motion Module D1
%*** * Cartesian inertial form, round rotating earth
%*** * Reserved C(3510) locations are 1600-1699
%*** * This module performs the following functions
%*** *
%*** * 1) Solves Newton's Law for spherical rotating earth in
%*** *    inertial coordinates
%*** * 2) Converts output to geographic variables
%*** *
%*** * MODIFICATION HISTORY
%*** * 960711 Created by Peter Zipfel
%*** *
%*** **************************************************************

  persistent accg agravg ai dum3 fspg fspv sbie sbii sbiid tei tge tgi tgv tig tvg vbeg vbei vbig vbii vbiid weii ; 
  grav=[];blon=[];blat=[];balt=[];dbi=[];dvbe=[];psivg=[];thtvg=[];dvbi=[];psivig=[];thtvig=[];
  global bcom_1; if isempty(bcom_1), bcom_1=zeros(1,3510); end;
  %% common c(3510);
  %% common bcom_1(3510);

  if isempty(fspg), fspg=zeros(1,3); end;
  if isempty(fspv), fspv=zeros(1,3); end;
  if isempty(agravg), agravg=zeros(1,3); end;
  if isempty(ai), ai=zeros(1,3); end;
  if isempty(tig), tig=zeros(3,3); end;
  if isempty(tei), tei=zeros(3,3); end;
  if isempty(sbie), sbie=zeros(1,3); end;
  if isempty(sbii), sbii=zeros(1,3); end;
  if isempty(vbei), vbei=zeros(1,3); end;
  if isempty(tge), tge=zeros(3,3); end;
  if isempty(tgi), tgi=zeros(3,3); end;
  if isempty(tvg), tvg=zeros(3,3); end;
  if isempty(tgv), tgv=zeros(3,3); end;
  if isempty(weii), weii=zeros(3,3); end;
  if isempty(vbig), vbig=zeros(1,3); end;
  if isempty(vbii), vbii=zeros(1,3); end;
  if isempty(vbiid), vbiid=zeros(1,3); end;
  if isempty(sbiid), sbiid=zeros(1,3); end;
  if isempty(dum3), dum3=zeros(1,3); end;
  if isempty(vbeg), vbeg=zeros(1,3); end;
  if isempty(accg), accg=zeros(1,3); end;

  %*** INPUT FROM EXECUTIVE

  % equivalence(bcom_1(0051),rearth) ::;
  % equivalence(bcom_1(0052),crad) ::;
  % equivalence(bcom_1(0058),weii3) ::;
  % equivalence(bcom_1(2000),t) ::;

  % CRAD = E Conversion from radians to degrees = 57.298
  % WEII3 = E Earth rotation - rad/sec

  %*** INITIALIZATION

  % equivalence(bcom_1(1622),tgv(1,1)) ::;
  % equivalence(bcom_1(1631),tig(1,1)) ::;
  % equivalence(bcom_1(1658),balt0) ::;

  % TGV(3,3) = I T.M. of  geographic wrt velocity coord - ND
  % TIG(3,3) = I T.M. of inertial wrt geographic coord - ND
  % BALT0 = I Saved value of initial altitude - m

  %***  INPUT FROM OTHER MODULES

  % equivalence(bcom_1(0205),grav) ::;
  % equivalence(bcom_1(1405),fspv(1)) ::;

  % GRAV= O Gravity acceleration - m/s^2
  % FSPV= O Specific force in velocity coordinates - m/s^2

  %*** STATE VARIABLES

  % equivalence(bcom_1(1640),vbiid(1)) ::;
  % equivalence(bcom_1(1643),vbii(1)) ::;
  % equivalence(bcom_1(1646),sbiid(1)) ::;
  % equivalence(bcom_1(1649),sbii(1)) ::;

  % VBIID(3) = S Time derivative of VBII(3) - m/s^2
  % VBII(3) = S Vel of missile wrt inertial frame in inertial axes - m
  % SBIID(3) = S Time derivative of SBIE(3) - m/s
  % SBII(3) = S Missile displacement from earth center in inertial axes - m

  %*** OUTPUT TO OTHER MODULES

  % equivalence(bcom_1(1606),balt) ::;
  % equivalence(bcom_1(1613),dvbe) ::;

  % BALT = O Vehicle altitude = m
  % DVBE = I/O Geographic speed - m/s

  %*** DIAGNOSTICS

  % equivalence(bcom_1(1602),psivgx) ::;
  % equivalence(bcom_1(1603),thtvgx) ::;
  % equivalence(bcom_1(1604),blon) ::;
  % equivalence(bcom_1(1605),blat) ::;
  % equivalence(bcom_1(1607),dvbi) ::;
  % equivalence(bcom_1(1608),psivigx) ::;
  % equivalence(bcom_1(1609),thtvigx) ::;
  % equivalence(bcom_1(1610),baltft) ::;
  % equivalence(bcom_1(1652),vbeg(1)) ::;
  % equivalence(bcom_1(1655),vbig(1)) ::;

  % PSIVGX = G Heading angle from north - deg
  % THTVGX = G Flight path angle from horizontal - deg
  % BLON = G Vehicle longitude - rad
  % BLAT = G Vehicle latitude - rad
  % DVBI = G Speed of vehicle wrt inertial frame
  % PSIVIGX = G Heading angle of inertial vel vect - deg
  % THTVIGX = G Flight path angle of inert vel vec  - deg
  % BALTFT = G Vehicle Altitude - ft
  % VBEG(3) = G Geographic velocity in geographic coord - m/s
  % VBIG(3) = G Inertial velocity in geographic coord - m/s

  %*** RIGHT HAND SIDE OF DYNAMIC EQUATIONS

  fspg*tgv,fspv,3,3,1);

  [agravg,dumvar2,dumvar3,grav]=vecvec(agravg,0.,0.,grav);

  [accg,fspg,agravg]=matadd(accg,fspg,agravg,3,1);

  [dum3,accg]=mateql(dum3,accg,3,1);

  ai*tig,dum3,3,3,1);

  %*** STATE VARIABLE INTEGRATION

  [vbiid,ai]=mateql(vbiid,ai,3,1);
  [sbiid,vbii]=mateql(sbiid,vbii,3,1);

  %*** UPDATE LONGITUDE, LATITUDE AND ALTITUDE, TVG AND FLIGHT PATH ANGLES

  [tei]=cadtei3(tei);
  sbie*tei,sbii,3,3,1);
  [blon,blat,balt,dbi,sbie]=cadsph3(blon,blat,balt,dbi,sbie);
  [tge,blon,blat]=cadtge3(tge,blon,blat);

  [weii]=matzer(weii,3,3);
  weii(1,2)=-weii3;
  weii(2,1)=weii3;
  dum3*weii,sbii,3,3,1);
  [vbei,vbii,dum3]=matsub(vbei,vbii,dum3,3,1);
  tgi*tge,tei,3,3,3);
  vbeg*tgi,vbei,3,3,1);
  [dvbe,psivg,thtvg,vbeg]=matpol(dvbe,psivg,thtvg,vbeg);
  psivgx=psivg.*crad;
  thtvgx=thtvg.*crad;

  %*** FOR NEXT INTEGRATION CYCLE: TIG, TGV

  [tig,tgi]=mattra(tig,tgi,3,3);
  [tvg,psivg,thtvg]=mat2tr(tvg,psivg,thtvg);
  [tgv,tvg]=mattra(tgv,tvg,3,3);

  %*** DIAGNOSTIC: INERTIAL VELOCITY IN GEOGRAPHIC AXES

  vbig*tgi,vbii,3,3,1);
  [dvbi,psivig,thtvig,vbig]=matpol(dvbi,psivig,thtvig,vbig);
  psivigx=psivig.*crad;
  thtvigx=thtvig.*crad;
  %---
  %-- Is this the right spot???
  %-- Save the altitude in feet.
  baltft=balt.*3.2808399;

  return;
end %subroutine d1